using Test
using FLEXINVERT
using NCDatasets
using Dates
using Printf

# Import the functions we want to test
using FLEXINVERT.Footprints: read_footprint_ncdf, find_release_window, save_averaged_footprint

@testset "Footprint Averaging Tests" begin

    @testset "find_release_window" begin
        # Test data: releases with start and end times in Julian days
        releases = zeros(Float64, 10, 2)
        base_jd = 2459000.0  # Some Julian day base

        # Create 10 releases, each 1 hour long, starting from base_jd
        for i in 1:10
            releases[i, 1] = base_jd + (i-1)/24.0  # start time
            releases[i, 2] = base_jd + i/24.0      # end time (1 hour later)
        end

        # Test case 1: Observation at start of releases, 4-hour averaging window
        obs_jdate = base_jd
        obs_avetime = 4.0/24.0  # 4 hours in Julian days
        start_idx, end_idx = find_release_window(obs_jdate, obs_avetime, releases)
        @test start_idx == 1
        @test end_idx == 4  # Should include releases 1,2,3,4

        # Test case 2: Observation in middle of release sequence
        obs_jdate = base_jd + 2.5/24.0  # 2.5 hours after start
        obs_avetime = 2.0/24.0  # 2 hour averaging window
        start_idx, end_idx = find_release_window(obs_jdate, obs_avetime, releases)
        @test start_idx == 4  # First release starting at or after obs time
        @test end_idx == 4    # Last release ending before obs + avetime

        # Test case 3: No releases in window
        obs_jdate = base_jd + 20.0/24.0  # 20 hours after start
        obs_avetime = 1.0/24.0
        start_idx, end_idx = find_release_window(obs_jdate, obs_avetime, releases)
        @test start_idx > end_idx  # No valid releases
    end

    @testset "NetCDF footprint I/O" begin
        # Create a temporary NetCDF file for testing
        tmpfile = tempname() * ".nc"

        # Test data dimensions
        nlon, nlat, ntime = 4, 3, 2

        # Create test coordinate data
        test_lon = Float64[0.0, 1.0, 2.0, 3.0]
        test_lat = Float64[50.0, 51.0, 52.0]
        test_time = Float64[2459000.0, 2459000.5]  # Two time steps

        # Create test grid data (time, lat, lon) as stored in NetCDF
        test_grid_raw = zeros(Float32, ntime, nlat, nlon)
        for t in 1:ntime, j in 1:nlat, i in 1:nlon
            test_grid_raw[t, j, i] = Float32(t * 100 + j * 10 + i)
        end

        # Write test NetCDF file
        NCDataset(tmpfile, "c") do ds
            defDim(ds, "longitude", nlon)
            defDim(ds, "latitude", nlat)
            defDim(ds, "time", ntime)

            lon_var = defVar(ds, "longitude", Float32, ("longitude",))
            lat_var = defVar(ds, "latitude", Float32, ("latitude",))
            time_var = defVar(ds, "time", Float64, ("time",))
            grid_var = defVar(ds, "grid", Float32, ("time", "latitude", "longitude"))

            lon_var.attrib["title"] = "longitude"
            lat_var.attrib["title"] = "latitude"
            time_var.attrib["title"] = "time"
            grid_var.attrib["title"] = "grid"
            grid_var.attrib["unit"] = "sm3/kg"

            lon_var[:] = test_lon
            lat_var[:] = test_lat
            time_var[:] = test_time
            grid_var[:, :, :] = test_grid_raw
        end

        # Test reading the file
        grid, gtime, lon, lat = read_footprint_ncdf(tmpfile)

        # Verify dimensions and data
        @test size(grid) == (nlon, nlat, ntime)  # Should be transposed to (lon, lat, time)
        @test length(gtime) == ntime
        @test length(lon) == nlon
        @test length(lat) == nlat

        # Verify coordinate data
        @test lon ≈ test_lon
        @test lat ≈ test_lat
        @test gtime ≈ test_time

        # Verify grid data is correctly transposed
        for t in 1:ntime, j in 1:nlat, i in 1:nlon
            expected_value = Float32(t * 100 + j * 10 + i)
            @test grid[i, j, t] ≈ expected_value
        end

        # Test round-trip: save and read back
        tmpfile2 = tempname() * ".nc"
        save_averaged_footprint(tmpfile2, Float32.(grid), gtime, lon, lat)

        grid2, gtime2, lon2, lat2 = read_footprint_ncdf(tmpfile2)

        @test grid2 ≈ grid
        @test gtime2 ≈ gtime
        @test lon2 ≈ lon
        @test lat2 ≈ lat

        # Clean up
        rm(tmpfile; force=true)
        rm(tmpfile2; force=true)
    end

    @testset "Footprint averaging integration" begin
        # Test that the main function handles configuration correctly
        using FLEXINVERT.Settings: Config
        using FLEXINVERT.ObservationsCore: Observations, ObservationRecord, add_observation!, empty_observations

        # Create test configuration with averaging disabled
        cfg = Config(average_fp=false)
        obs = empty_observations()
        files = nothing  # Won't be used since averaging is disabled

        # This should return early without error (info log only)
        @test_logs (:info, r"average_fp disabled") FLEXINVERT.Footprints.average_footprints!(obs, cfg, files)

        # Test with averaging enabled but empty observations
        cfg_enabled = Config(average_fp=true)
        @test_logs (:info, r"Starting footprint averaging") (:info, r"Footprint averaging completed") begin
            FLEXINVERT.Footprints.average_footprints!(obs, cfg_enabled, files)
        end
    end

end
